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1.  OVERVIEW  AMD  SUMMARY 


Fundamental  to  tfaia  work  ia  the  development  of  a  continuum 
formulation  that  can  accurately  account  for  the  effects  oi 
interlaminar  shear  and  interlaminar  normal  stress  variation 
thru-the-thickness  of  a  laminate.  Furthermore,  emphasis  is 
particularly  on  tapered-twisted  airfoil  geometries  which  can  be 
analytically  represented  as  an  assemblage  of  thin  to  moderately  thick 
finite  elements.  To  achieve  solution  efficiencies,  the  elements 
developed  in  this  work  are  of  the  triangular/ quadrilateral  plate  type 
as  opposed  to  solid  type  elements. 

On  the  basis  of  these  requirements  and  considering  viable 
alternatives,  three  suitable  continuum  formulations  have  been 
developed  and  are  herein  denoted  as  the  (i)  Higher  Order  Displacement, 
(ii)  Modif ied-Kirchhot t  and  (iii)  Hybrid  Stress  formulations, 
respectively.  The  former  two  formulations  have  been  incorporated  in  a 
computer  code  and  the  various  elements  have  been  tested  on  the  basis 
of  correlations  with  known  analytical,  numerical,  and  experimental 
solutions.  Numerous  tests  have  been  performed  for  linear  static  and 
linear  dynamic  cases.  It  is  noted  that  the  code  has  some  unique 
features,  e.g.,  it  can  assemble  elements  having  an  unequal  number  of 
degrees  of  freedom  at  its  nodes,  it  treats  arbitrary  ply  orientations 
and  it  performs  integration  on  a  layer-by-layer  basis  through  the 
laminate.  Herein  a  layer  refers  to  either  a  lamina  or  to  a  sub-set  of 
laminae  having  equal  ply  orientations.  The  latter  feature  is 
essential  in  developing  a  fully  nonlinear  capability. 


Significant  efforts  have  also  been  devoted  to  neveloping  a 
suitable  large  displacement  formulation.  Due  to  the  requirement  that 
interlaminar  stresses  be  accurately  represented ,  a  total  Lagrangian 
formulation  is  utilized  and  is  based  upon  the  complete  Green's  strain 
tensor.  A  geometric  and  large-aisplacement  stiffness  formulation  has 
been  implemented  in  the  computer  cone  baaed  upon  a  form  of  the 
nonlinear  strain-nodal  displacement  relationships  suitable  for  each  of 
the  elements  under  development. 

An  extensive  literature  survey  has  been  performed  to  identify 
analytically  tractable  methods  of  treating  damage  accumulation  in 
composites.  Since  emphasis  in  this  work  is  on  the  development  of 
incremental  response  solutions,  the  computational  approach  must  have 
the  capability  to  (i)  predict  and  differentiate  between  relevant 
failure  modes,  (ii)  modify  constitutive  equations  appropriately  and 
(iii)  perform  equilibrium  iterations  to  assure  stress  redistribution 
based  upon  the  extent  of  damage.  Use  of  "piecewise  smooth"  failure 
criteria  based  on  various  types  of  damage  provides  a  good  basis  for 

e 

incrementally  tracking  damage.  Ibis  approach  has  been  incorporated  in 
the  computer  code.  Mote  that  integration  for  an  element  is  performed 
on  a  layer-by-layer  baais  which  allows  tor  damage  effects  to  be 
characterized  at  the  layer  level.  It  is  noteworthy  that  variation  in 
strain  energy  can  be  calculated  as  damage  accumulates  and  that  it  may 
be  possible  to  go  further  to  predict  useful  strain  energy  release  rate 
values.  Ihus  it  may  be  possible  to  make  use  of  energy  in  addition  to 
maximum  stress  criteria  to  characterize  damage. 

Experimental  data  of  the  type  required  to  substantiate  damage 
predictions  has  been  assembled  to  the  extent  possible.  Analysis/test 


correlations  have  been  performed  for  selected  laminates.  It  is  notea 
that  useful  experimental  data  is  quite  limited. 

Technical  progress  in  this  program  has  been  subs  tan"' -illy  on 
schedule  vith  regard  to  developing  continuum  formulations  and  a 
suitable  finite  element  code.  It  has  not  been  possible,  however,  to 
complete  the  damage  characterization  efforts.  Work  will  continue 
under  an  extension  to  this  contract  with  the  intent  of  fully 
implementing  damage  characterization  in  the  nonlinear  transient 
analysis. 


considerations.  Finite  element  model*  have  been  teetea  tor  linear 
static,  dynamic  and  buckling  analysis,  lhe  teat  problems  and  tbe 
results  are  presented  in  Section  11.1.4.  The  finite  element  models 
are  herein  briefly  discussed. 

A.  Higher. Order  Displacement  Formulation 


lhe  thru- the- thickness  effects  can  be  incorporated  into  an 
analysis  by  choosing  a  displacement  field  that  eliminates  two  major 
shortcomings  of  the  classical  plate  theory;  namely  normals  remain 
normal  and  in-plane  displacements  are  linear  thru  the  thickness, 
these  shortcomings  are  eliminated  by  prescribing  indepenaently  the 
reference  surface  displacements  ana  rotations  of  the  normal  and 
including  higher  oraer  terms  for  in-plane  displacements,  this  iB 
accomplished  by  the  following  variation 

u(x»yiZ)  -  u0(x,y)  +  Z-X(x,y)  +  z2^x(x,y) 
v(x,y,z)  =  v0(x,y)  +  z;y(x,y)  +  z 2iy(x,y) 
w(x,y,z)  =  w0(x,y) 


tbe  neutral  surface  aisplacements  are  represented  by  u  ,  v  and  w  , 

0  0  0 

the  rotation  about  y-axiB  is  denoted  by  ip*  and  the  rotation  about  the 

z 

x-axis  is  .  the  coefficients  of  z  ,  i.eM  <p  and  4>  ,  are 

y  x  y 

contributions  from  transverse  deformations  [5,b] . 


tbe  elements  developed  are  designated  as  the  quadrilateral  higher 
order  displacement  (QhD)  models.  QHD40  is  an  eigbt-noded  element  with 
seven  degrees  of  freedom  (three  midsurface  displacements,  two 
rotations  and  two  higher  order  terms  for  in-plane  displacements)  per 
corner  node  and  three  degrees  of  freedom  (transverse  midsurface 
displacement  and  two  rotations)  per  mid-side  node.  Element  QHD28  is 


6 


simplified  version  of  QHD40  vbere  tbe  mid-side  nodes  are  eliminated 


/ 


element  with  27  d.o.f.-  TD27M.  The  stress  fields  obtained  for  these 
elements  represents  a  quadratic  thru  the  thickness  variation  for  the 
transverse  shear  stresses  and  a  cubic  variation  for  the  transverse 
normal  stress.  The  respective  displacement  fields,  noda1  degrees  of 
freedom  and  stress  fields  are  given  in  Appendix  IB. 

11.1.2.  Large  Displacement  Formulation 

Inclusion  of  geometrically  nonlinear  effects  in  the  formulation 
must  be  based  upon  both  the  geometry  to  be  analyzed  and  upon  the  type 
of  stress  prediction  capabilities  desired.  The  classical  approach  to 
thin  plate  analysis  has  been  to  use  the  Kirchhof f-Love  assumptions  in 
conjunction  with  the  nonlinear  von  Harman  relations  [ll,12j.  As 
previously  indicated,  the  Kirchhof t-Love  assumptions  are  relaxed  in 
this  vork  to  allow  tor  a  more  accurate  definition  of 
interlaminar-shear  and  interlaminar-normal  stress  variations.  These 
stresses  can  vary  substantially  through-the-tbickness  for  tbe 
geometries  of  interest,  i.e.,  thin  to  moderately  thick  plate  type 

4 

structures.  Furthermore,  the  requirement  that  these  stresses  be 
accurately  determined  means  that  tbe  nonlinear  portion  of  tbe 
strain-displacement  relationship  must  contain  all  significant 
coordinate  displacements.  The  complete  Green's  strain  tensor  is 
utilized  in  this  work,  therefore,  to  account  for  all  significant 
contributions  to  the  interlaminar  stress  field,  kith  respect  to 
fixed  Cartesian  coordinates,  x,  y,  and  z,  tbe  strain  tensor  has  the 


♦  ' 


form 


8 


H  +  I  f/Ui)2  +  (ll)2  +  /3w^l 

3x  2  l'3x'  l3x'  'ax'  j 


Y  *  iii  +  +  fiii  l!i  +  lx  3v_  .  3w  aw 

Txy  3y  3x  1 3x  3y  ax  3y  37  ay 


where  u,  v  and  w  represent  displacements  in  the  x,y,a  coordinate 


directions,  respectively*  Note  that  the  other  strain  components  are 


obtained  by  a  suitable  permutation.  In  small-displacement  analysis. 


the  quadratic  terms  are  neglected  to  give  simply  the  linear  strain 


approximation. 


Based  on  the  Greek's  strain  tensor,  the  strain  to  nodal  point 


displacement  relationship  can  be  specified  for  elements  under 


development.  It  takes  the  form 


{e}  =  (B]{d} 


where  is  the  vector  of  strain  components,  {£}  the  vector  of 


nodal  point  displacements  and  [B]  a  function  of  derivatives  of  the 


element  shape  functions.  The  quadratic  terms  in  the  strain  tensor 


result  in  IB]  being  a  function  of  displacement  state  and,  therefore, 


an  incremental  equilibrium  formulation  is  required.  Ibe  incremental 


strain-nodal  displacement  relationship  takes  the  form 


Ue}  =  ([Bq]  +  [Bl])  C 5 a } 


where  {<$e|  and  {oA}  represent  incremental  strains  and  nodal 


displacements,  respectively,  [B  J  and  [B  ]  are  the  small  and  large 

O  L 


displacement  contributions  to  the  incremental  strains.  Based  on  the 


incremental  equilibrium  equations,  the  displacement  formulation  gives 
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the  force-displacement  relationships 
(K01  -  /  (S„]T  [OHB0]  <JV 

V 

tKL3  “  /  K^T  [DJt8L]  +  C8l3T[D]  [Bl]  +  [Bl]T[D]  [B0j  y  dV 

v 

where  IDJ  is  an  elasticity  matrix  obtained  simply  from  the 
constitutive  equations  and  integration  is  over  the  volume  V  of  the 
element.  IKq1  is  denoted  the  small-displacement  stiffness  matrix  and 
tK^]  is  denoted  the  large-displacement  stiffness  matrix.  Since 
response  is  also  a  function  of  stress  state,  the  geometrical  stiffness 
matrix  [Kq]  is  required  and  is  obtained  from 

fKGJ(6A}  =  /-S(8LIT{0|dV 

v 

where  jcj  is  the  vector  of  stress  components. 

Inertial  effects  are  analytically  treated  as  a  mass  matrix  1M] 
which  is  a  function  of  density  and  the  element  shape  functions  (see 
Appendix  II).  Ihese  matrix  forms  are  required  in  formulating 

static/dynamic  response  solutions  and  the  incremental  equilibrium 
equations  have  the  general  fo:.*m 

Cn]t6u>  +  ([K01  +  [kl3  +  Ckg])  (*u>  = 

where  the  use  and  stiffness  matrices  represent  an  assembly  of  the 
elemental  matrices  previously  discussed,  |$u}  and  {<5u|  represent  the 
incremental  displacements  and  accelerations  tor  the  mathematical  model 
and  |5f|  represents  the  vector  of  incrementally  applied  forces. 

In  developing  a  geometrically  nonlinear  formulation,  the  effort 
is  largely  in  defining  the  incremental  strain-nodal  displacement 


relationship.  Having  developed  this  relationship  for  a  particular 
element,  stiffness  matrices  are  readily  developed  as  the  preceding 
equations  indicate.  These  relstionships  are  presented  :  'ppendix 
III.  The  form  of  these  equations  is  the  same  for  all  elements. 

II. 1.3.  Computer  Implementation 

A  computer  code  has  been  developed  for  the  purpose  of 
implementing  the  various  continuum  formulations.  At  present ,  the  code 
performs  the  following  fundamental  calculations: 

element  stiffness  matrix  (linear,  nonlinear,  geometric) 
generation 

element  mass  matrix  generation 
assembly  of  equilibrium  equations 
decomposition  and  solution  of  equilibrium  equations 
equilibrium  iteration  for  incremental  solutions 
fundamental  frequency  and  mode  shape  calculation 
elastic  buckling  calculation 

A  characteristic  of  the  elements  under  development  is  that  node 
points  can  have  different  numbers  of  degrees  of  freedom,  i.e., 
typically  mid-side  nodes  have  fewer  degrees  of  freedom  than  corner 
nodes.  The  code  has  been  fashioned  to  handle  this  condition.  All  of 
the  integration  is  performed  on  a  layer-by-layer  basis  thru  the 
thickness  of  the  laminate.  This  approach  is  fundamental  to  developing 
the  capability  to  allow  for  inelastic  material  behavior  and, 
ultimately,  to  the  inclusion  of  damage  mechanisms  in  the  formulation. 


Since  solution  of  tbe  equilibrium  equetiona  ia  a  vital  component 

in  tbe  overall  aolution  strategy,  it  ia  appropriate  to  discuss  the 

numerical  methodology  uaed  in  aolving  tbeae  equations.  Ibe  intent  is 

to  obtain  a  higher  ordered  variation  of  the  transverse  shear  and 

normal  stresses  (a  ,  a  and  a  )  than  can  be  obtained  v  ae 
xz  yz,  zz 

equilibrium  equations.  The  solution  procedure  can  be  thought  of  as 

described  belov.  Assume  that  the  in-plane  stresses  (<j  <j  a  ) 

xx  yy  ‘  xy 

within  each  layer  of  a  particular  element  have  been  determined  at 
selected  locations,  i.e.,  through  solution  of  the  constitutive 
equations.  In  the  code  as  presently  written,  these  locations  are 
specified  as  tbe  element  centroid  and  element  nodal  points.  The 
equilibrium  equations  (in  the  absence  of  body  forces)  have  the 
indicial  form 

*  0 

from  which  it  follows  that  the  tbru-the-thickness  shear  stress 
variation  can  be  written  in  numerical  form  tor  the  layer  as 

lc*Zi  =  -( 'xx.x  +  '-rxy»y/-j  M 

and 

L°W  =  ’(jxy,x  +  -yy.y),  tZi 

Here,  the  left-band-aide  represents  the  change  in  stress  from  the 
lower  to  the  upper  surface  of  tbe  1thlayer  and  AZ£  is  the  thickness  of 
the  ithlayer  at  a  particular  location.  The  derivatives  with  respect 
to  x  and  y  in  tbe  expressions  above  are  readily  computed;  this  is 
because  in-plane  stresses  within  a  layer  are  related  to  element 
displacements  through  derivatives  of  element  shape  functions  in 
conjunction  with  a  material  aefinition. 

For  an  n  layered  laminate,  n  equations  can  be  written  in  terms  of 


both  the  unknown  sheet  stresses  et  layer  interfaces  and  the  shear 
stresses  at  the  laainate  surfaces.  Assuming  the  laainate  has 
shear-free  surfaces,  the  equations  above  give  n  equations  in  n-1 
unknowns,  so  that,  the  equation  set  is  over-determined.  .  equations 
have  eke  matrix  form  below 


"  x  <n-l)  (n-1)  x  1  (n  x  i) 


where  IX2i3  *  (0xx,x  +  0xy,y)^  and  0^  represents  the  shear  stress 
acting  at  the  interface  of  the  j-1 th  and  jth  layer.  A  similar 
equation  set  is  obtained  by  replacing  cx2j  with  cyz;j  and  Ixzi  with  I  t 
These  equations  are  solved  by  utilizing  a  least-squares 
orthonormalization  procedure  [133.  Due  to  the ’simplicity  of  the  terms 
in  the  coefficient  matrix,  a  concise  closed-form  solution  is  obtained. 
Having  determined  the  transverse  shear  stresses,  the  transverse 
normal  stress  variation  is  determined  through  the  numerical  form  of 
the  third  equilibrium  equation  for  the  i th  layer 

“JZZ}  =  -(axz’X  +  ayz»y)i  Z.Z. 

As  before,  the  left-hand-side  represents  the  change  in  stress  through 
the  1th  layer.  Appropriate  polynomial  functions  are  utilized  to 
describe  the  a  and  ov_  in-plane  variation.  These  functions  are 


differentiated  to  obtain  the  right-hand-side  of  the  equations  above. 
Again  the  equation  set  is  overdetermined  because  the  normal  tractions 
are  known  at  the  laminate  surfaces.  Solving  for ozz  proceeds, 
therefore.  In  identically  the  same  manner  as  discussed  in  calculating 
axz  and  oyz.  Parenthetically,  inclusion  of  body  forces  a.  -  later  date 
can  be  accomplished  with  little  difficulty. 

It  should  be  emphasized  that,  the  successful  application  of 
Higher  Oroer  Displacement  type  elements,  i.e.,  for  particularly  thin 
geometries  is  to  utilize  reduced  numerical  integration  where  as  this 
is  not  necessary  for  the  Modified  Kircbhoff  formulation.  This 
approximation  technique  brings  along  the  choice  of  implementing  it 
overall  or  selectively  to  the  strain  energy  components.  For  the  QHD 
formulation,  only  the  transverse  shear  components  are  integrateo  with 
reduced  order  [14-16].  An  undesirable  aspect  of  this  approach  is  that 
the  reduced  integration  order  may  affect  the  physical  behavior  of  the 
element  by  introducing  spurious  zero  energy  modes.  It  is  desirable  to 
have  only  rigid  body  modes  since  there  does  not  yet  seas  to  be  a 
generally  accepted  method  of  controlling  the  additional  modes. 

II. 1.4.  Analytical  Verification 

As  noted,  elements  formulated  on  the  basis  of  independent 
transverse  displacements  and  rotations,  require  reduced  quadrature  for 
good  performance.  For  QHD4Q,  3x3  Gaussian  quadrature  along  with  the 
2x2  quadrature  for  the  transverse  shear  components  is  employed.  QHD28 
and  QHD20  formulations  are  similarly  integrated  with  2x2  ano  lxl 
Gaussian  quadratures.  Manipulation  of  quadrature  rules  may  produce 
spurious  zero  energy  modes  in  addition  to  the  required  rigid  body 


■odea,  thus  detracting  from  overall  element  performance.  116,17].  A 
apectral  (eigenvalue)  test  has  been  conducted  with  and  without  full 
quadrature  to  observe  the  zero  energy  modes  of  the  QHD  ments.  The 
quadrature  order,  the  number  of  zero  eigenvalues  and  tn«  corresponding 
number  of  spurious  zero  energy  modes  for  the  QUD4U,  QRD28 ,  and  QHD20 
elements  are  listed  in  Table  1.  The  spurious  mode  shapes  associated 
with  the  QHD28  element  are  illustrated  in  Figure  1.  Since  the  QHD40 
formulation  does  not  exhibit  spurious  modes,  it  can  be  utilized  in 
modelling  complex  geometries  without  concern  for  controlling  such 
behavior. 

It  is  also  noteworthy  to  observe  the  effect  of  reduced 
integration  on  the  representation  of  the  generalized  forces.  In 
order  to  illustrate  the  effect,  the  forces  associated  with  the 
transverse  displacement  of  a  corner  node  are  sketched  in  Fig.  2  tor 
QU026  with  and  without  reduced  integration  respectively. 

In  the  examples  that  follow,  performance  of  the  QHD  formulation 
is  demonstrated  by  comparing  results  to  those  obtained  by  classical 
plate  theory  (CFT),  by  elasticity  and  by  other  finite  element 
formulations  for  linear  static  dynamic  and  buckling  analyses. 

Limited  results  are  also  presented  for  the  QD  and  ID  formulations. 

The  latter  results  are  simply  presented  for  comparison  because  it  is 
apparent  that  the  QhD  formulation  always  gives  the  best  results.  It 
would  seem,  therefore,  that  the  higher  order  displacement  formulation 
is  the  better  approach.  The  ortho  tropic  material  properties  used 
throughout  are  tabulated  in  Table  2.  Geometries  studied  include 
cylindrical  bending  of  a  plate  as  well  as  bending  of  simply  supported 


square  and  rectangular  plates.  Various  ply  layups  are  considered  and 
loading  is  that  of  a  sinusoidally  and  uniformly  distributed  pressure. 
Cylindrical  bending  is  modelled  as  a  strip  of  twenty  elements.  For 
the  simply  supported  plates,  symmetry  considerations  allow  that  only  a 
quadrant  of  the  plate  need  be  modelled.  Fineness  of  th  .b  is 
varied  to  demonstrate  solution  convergence.  Additionally,  distorted 
meshes  are  considered  to  demonstrate  modelling  considerations.  For 
the  examples  involving  symmetric  layups,  the  quadratic  terms  of  QHD4Q 
are  restrained;  so  that,  32  degree  of  freedom  elements  are  utilized  to 
obtain  these  solutions.  This  is  allowable  in  these  particular  cases 
because  the  quadratic  terms  do  not  significantly  affect  the  results. 
This  is  not  true  in  the  first  example  considered. 

Static  Response  Calculations 

Cylindrical  Bending-  Bidirectional  (0  / 90)  Sine  Load,  Material  11 

Fibers  run  parallel  to  the  plane  of  curvature  in  the  lower  layer 
and  are  rotated  90°  in  the  upper  layer  of  the  plate.  Layers  are  of 
equal  thickness  which  is  true  in  the  subsequent  example  problems  as 
well.  Ihe  elasticity  solution  obtained  by  Pagano  [9J  gives  a  nearly 
quadratic  z  variation  in  u,  where  u  is  the  normalized  in-plane 
displacement  of  the  laminate.  In  this  instance,  inclusion  of  the  z2 
terms  in  the  finite  element  modelling  should  affect  the  results.  This 
is  demonstrated  in  Figures  3  to  5.  Results  demonstrate  differences 
obtained  with  and  without  quadratic  terms.  The  difference  is  greatest 
for  the  lower  aspect  ratios,  e.g.,  tor  S  ■  4  a  difference  of  111  is 
obtained.  In  Figure  4,  the  calculated  normalized  in-plane  stress 
variation  is  presented  tor  an  aspect  ratio  of  4.  Note  that  maximum 
variation  is  presented  for  an  aspect  ratio  of  4.  Note  that  maximum 
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stresses  differ  by  some  36X  when  computed  with  sad  without  the  z 
terms,  respectively*  The  effect  of  including  quadratic  terms  in  the 
finite  element  solution  is,  therefore,  much  more  pronoun' ad  when 
stresses  as  opposed  to  displacements  are  considered.  i-_~;e  5 
demonstrates  this  effect  on  stress  computation  as  a  function  of  aspec 
ratio.  Mote  that  calculated  quantities  are  normalized  in  this  example 
and  in  those  that  follow  as  in  the  cited  references. 


Cylindrical  Bending-Symmetric  (0  / 90  /0)  Sine  Load,  Material  II 

For  this  geometry,  fibers  are  parallel  to  the  plane  of  curvature 
in  the  outer  layers  and  rotated  90°  in  the  middle  layer.  Calculated 
stresses  are  compared  to  the  elasticity  solution  of  Pagano  [9]. 
Figures  6  and  7  present  the  normalized  transverse  shear  stress 
variation  at  the  6imply  supported  boundary.  Figures  8  and  9  present 
the  normalized  in-plane  stress  variation  at  the  center  of  the  bent 
surface. 


Simply  Supported  Square  Plate  (0  / 90  /0)  Sine  Load,  Material  11 

Fibers  in  the  outer  layers  of  the  laminate  run  parallel  to  the  x 
axis  while  those  in  the  middle  layer  run  parallel  to  the  y  axis,  where 
the  origin  of  coordinates  is  located  at  a  corner  of  the  plate  and  in 
the  mid-plane  (see  Figure  107.  This  coordinate  system  is  consistent 
with  examples  that  follow  as  well.  Consider  the  plate  as  having 
planar  dimensions  a  x  a  and  total  thickness  h.  Solutions  have  been 
generated  for  aspect  ratios  S  ■  4  to  100,  where  S  ■  a/h.  Transverse 
shear  stress  variation  <?xz  at  (x,y)  coordinates  (0,a/2)  and  in-plane 
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shear  stress  variation  a  at  coordinates  (0,0)  are  presented  in 

xy 

Figures  10  and  11  tor  an  aspect  ratio  of  4.  Mote  that  the  comparison 

is  between  the  present  finite  element  results  and  those  obtained  via 

elasticity  (l&l  and  CFT.  Calculated  short-transverse  nc  stress 

variation  o  ii  presented  in  Figures  12  and  13.  These  stresses  are 
z  z 

normalised  as  a__«a,/10Q  at  the  center  of  the  plate  and  as  a  ■  10a 

ZZ  ZZ  2.2,  ZZ 

at  the  edge  of  the  plate.  Results  are  compared  to  those  obtained  by 
elasticity  over  a  range  of  aspect  ratios  in  Table  3.  Similar  results 
are  given  in  Table  3.1  tor  the  QD  formulation  and  in  Table  3.2  tor  the 
TD  formulation  .  Convergence  characteristics  are  demonstrated  by 
presenting  results  obtained  using  2x2,  3x3,  and  6x6  meshes.  The  finer 
mesh  gives  better  agreement,  but  the  coarser  mesh  gives  very 
reasonable  correlation  also. 

The  effects  of  distorting  the  mesh  have  also  beeu  considered  to  a 
limited  extent.  Results  have  been  obtained  for  the  relatively  coarse 
meshes  shown  in  Figure  14.  Calculated  stresses  and  displacements  are 

a 

presented  as  a  function  of  aspect  ratio  and  compared  to  the  elasticity 
solutins  in  Table  4.  As  expected,  the  values  are  not  as  accurately 
determined  as  are  those  obtained  via  the  regular  meshes.  Distortion 
of  the  mesh  has  a  much  more  dramatic  effect  upon  the  calculated 
transverse  shear  stresses  than  upon  the  calculated  in-plane  stresses 
and  displacements.  Since  the  transverse  stresses  are  based  on 
equilibrium  considerations,  it  seems  the  mesh  must  be  refined  enough 
to  reasonably  approximate  equilibrium.  This  is  especially  apparent  in 
comparing  results  obtained  tor  mesh  A  to  those  obtainea  for  mesh  C. 

1°  each  of  these  cases,  elements  having  a  taper  ratio  of  2  to  1  are 
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utilized.  Mesh  C  gives  significantly  improved  transverse  stresses, 
however,  because  the  mesh  is  tine  enough  to  better  represent  the 
loading  distribution. 

Simply  Supported  Rectangular  PlateCO  / 90  /0)  Sine  Loao,  Material  11 
Ortbotropic  layers  have  the  same  orientation  as  in  the  previous 
example.  The  plate  has  dimensions  a  x  b,  where  b  is  three  times  a. 
Solutions  have  been  obtained  for  aspect  ratios  (S  ■  a/h)  ranging  from 
4  to  100.  Transverse  shear  stress  variation  Qyz  coordinates  (a/2,0) 
is  given  in  Figure  15  tor  an  aspect  ratio  S  ■  4.  Comparison  is  made 
to  both  elasticity  and  CPT  solutions.  A  full  range  of  results  are 
presented  in  Table  5  and  compared  to  those  obtained  via  elasticity 
118]  and  to  those  obtained  by  Reddy  [IV]  in  a  recent  finite  element 
formulation.  Correlation  with  elasticity  is  quite  good,  particularly 
tor  aspect  ratios  of  lb  and  above,  and  appear  to  be  more  accurate  than 
those  obtained  with  the  alternate  finite  element  solution.  Results 
obtained  using  the  QD  and  TD  elements  are  poor  (not  presented  herein) 
compared  to  those  obtained  using  the  QHD  element.  Thus  the  QD  and  TD 
elements  are  quite  sensitive  to  element  distortion  and  not  suitable 
for  general  analysis  purposes. 

Simply  Supported  Square  Plate  (0  / 90  /V0  /0)  Sine  Load,  Material  11 
The  laminate  geometry  consists  of  outer  layers  with  fibers 
parallel  to  the  x  axis  and  inner  layers  with  fibers  parallel  to  the  y 
axis.  The  plate  has  planar  dimension  a  x  a  and  total  thickness  b. 
Stress  and  displacement  results  are  presented  in  Table  6  for  aspect 
ratios  ranging  from  4  to  100.  Similar  results  are  given  in  Table  6.1 


based  on  the  QD  formulation  and  in  Table  6.2  based  on  the  lb 
formulation.  Results  are  compared  to  both  elasticity  and  to  other 
finite  element  results.  Again,  the  computed  values  are  in  excellent 
agreement  with  elasticity  120]  for  moderately  thick  to  'n  geometries 
and  are  more  accurate  than,  the  compared  to  numerical  x  ...its. 

Solutions  also  have  been  obtained  for  the  present  geometry  on  th. 
basis  of  reduced  vs.  full  integration.  This  comparison  is 
demonstrated  in  Figures  16  and  17  by  giving  percent  error  in 
calculated  values  vs.  aspect  ratio.  It  is  apparent  that  reduced 
integration  is  particularly  needed  to  minimize  errors  in  calculated 
transverse  stresses  and,  furthermore,  solution  validity  over  a  vide 
range  of  laminate  geometries  is  demonstrated. 

Fundamental  Frequency  Calculations 

To  assess  the  effects  of  finite  element  formulations,  aspect 
ratio,  support  conditions  and  the  lamina  stacking  sequences  on  the 
fundamental  natural  frequencies  of  composite  plates,  the  problems 
listed  in  Table  7  are  considered  [21] . 

The  nonrdimensionalized  fundamental  frequency  tor  the  cross-ply 
laminate  of  Problem  1  versus  aspect  ratios  is  given  in  Table  8.  As 
can  be  seen,  all  three  elements  predict  frequencies  that  are  in 
excellent  agreement  with  the  closed  form  solutions  obtained  by  Reddy 
122]. 

The  effects  of  higher  order  terms  in  the  displacement  based 
finite  element  formulations  are  investigated  for  Problem  2.  here,  the 
performances  of  QHD40  and  QHU28  (with  higher  order  terms  locked)  are 
compared  to  elements  STPD1  and  STPD3  of  123]  with  linear  and  cubic 
v*r*'*tions  through  the  thickness  respectively.  The  results  are 


summarized  in  Table  9.  The  normalized  fundamental  frequencies  of 
Problem  3  are  displayed  in  Figure  18.  Mote  that  the 
non-dimenaionalized  fundamental  frequency  increases  as  the  angle  of 
orientation  is  increased  for  both  symmetric  and  antisy  ric 
angle-ply  square  plates.  This  observation  is  in  excellent  agreement 
with  Reddy/ s  [22]  antisymmetric  laminate.  In  Figure  19,  a  decrease  in 

N. 

the  fundamental  frequency  is  observed  as  the  angle  of  orientation  is 
increased  for  the  angle-ply,  cantilever,  rectangular  and  square  plates 
of  Problem  4.  The  difference  between  Figures  18  and  19  are  attributed 
primarily  to  the  different  support  conditions. 

Further  investigations.  Problem  3,  of  angle-ply  laminates  are 
summarized  in  Table  10.  The  stacking  sequences  of  reference  [24J  are 
used  to  illustrate  their  effects  on  the  fundamental  frequency 
calculations.  The  numbers  within  parenthesis  are  calculated  by 
Crawley  124,  25]. 

Shown  in  Figure  20  is  the  variation  of  the  non-dimenaionalized 
fundamental  frequency  for  cylindrical  bending  problem,  as  calculated 
via  the  QD  formulation  and  the  classical  plate  theory.  For  comparison 
purposes,  the  frequencies  are  normalized  with  respect  to  the  classical 
plate  theory  results. 

Transient  Response  Calculations 

Element  performance  has  been  evaluated  with  respect  to  predicting 
linear-transient  response.  Both  displacements  and  stresses  have  been 
determined  tor  a  variety  of  laminated  plate  geometries  subjected  to 
instantaneously  applied  pressure  loading.  These  results  have  been 
compared  to  those  obtained  via  both  CPI  and  a  shear  deformable  theory 
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(SD1)  [2b]  .  Typical  results  are  presented  in  Figures  21  ana  22  tor  a 
(0/90)  square  plate.  In  this  example,  the  plate  is  quite  thick  in 
that  it  has  an  aspect  ratio  ot  3. 

Buckling  Calculations 

Element  pertormance  has  also  been  veritiea  tor  linear  elastic 
buckling  calculations.  As  an  example,  critical  buckling  loan  is 
plotted  vs.  the  E^/E2  ratio  tor  anti-symmetric  (0/90)  and 
(0/90/0/90/0/90)  cross  ply  composite  plates  in  Figure  23.  kesults 
agree  well  with  those  obtained  by  boor  127 J  and  with  other  shear 
detonuable  theories.  Good  results  have  also  been  obtained  tor  other 
boundary  conditions  and  stacking  sequences. 
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11.2.  TASK  11:  Incorporate  Damage  Mechanisms  into  Dynamic  Response 
Formulation 

The  literature  survey  1 2b— 63 J  pertormeo  has  been  quite  helptul  in 

terms  ot  delineating  the  viable  approaches  to  including  damage 
mechanisms  in  the  analysis.  Relevant  tailure  modes  ct  interest 

include  those  listed  below 

(l)  tiber  tracture 

(li)  liber-matrix  debonaing 

(iii)  matrix  cracking  (parallel  and  transverse  to  tibers) 

(iv)  oelamination 

(v)  buckling  (possibly  at  layer  or  sub-laminate  level) 

Several  smooth  tailure  criteria,  e.g.,  L 64—67 J  have  been  developed  in 
recent  years  to  represent  the  tailure  ot  composites.  These  criteria, 
to  varying  degrees,  can  predict  "tailure"  but  uo  not  identity  a 
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particular  mode  of  failure.  1b  performing  incremental  "damage" 
analyais,  it  ia  essential  to  both  predict  failure  and  to  characterize 
itt  e.g.,  do  fibers  rupture,  does  delamination  occur*  etc.  The 
computational  approach  must*  therefore,  differentiate  between  viable 
failure  modes  and  appropriately  alter  the  constitutive  equations  on 
an  incremental  basis.  This  can  be  accomplished  by  impementing  a 
piecewise  smooth  failure  criteria,  e.g.,  128]  in  the  finite  element 
formulation.  The  general  failure  criteria  is  then  comprised  of  m 
separate  inequalities  of  the  form 

Fj  ({<?})«  1  ;  j  =  1,2,. ...m 

at  the  layer  level  within  each  element.  These  criteria  should 
differentiate  between  (i)  tensile  and  compressive  fiber  failure,  (ii) 
tensile  and  compressive  matrix  failure  and  (iii)  oelamination  at 
layer  interfaces  due  to  either  maximum  stress  or  buckling 
considerations. 

As  progressive  damage  occurs  throughout  incremental  loading 
(whether  it  be  static  or  dynamic),  it  is  essential  that  violation  of 
failure  criteria  inequalities  be  reflected  in  modification  of  the 
material  properties.  This  can  be  achieved  by  including  damage  state 
variables  1 47]  in  the  constitutive  equations  to  reflect  "stiffness 
reduction."  These  equations  can  be  represented  as 

io}  =  [D][y]u} 

where  ID]  represents  the  material  matrix  and  [y  ]  contains  the  damage 
state  variables.  The  latter  provide  the  basis  for  changing  the  O^j 
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tens  based  upon  tbe  extent  to  which  the  failure  criteria  are 
violated. 

In  conjunction  with  the  above  it  is  essential  to  perform 
equilibrium  iterations  within  each  analysis  increment.  Ibis  is 
required  to  assure  that  stress  redistribution  is  properly  accounted 
for  as  damage  progresses. 

Tbe  failure  criteria  currently  implemented  in  the  incremental 
analysis  are  primarily  due  to  Hashin  [28],  Lee  129J,  Greszczuk  [31] 
and  Habn  139].  Layer  stresses  are  defined  as  shown  in  Figure  24  and 
the  criteria  are  summarized  as  follows: 


Fiber  Failure 


1 .  Tension 

The  simplest  criterion  for  tensile  failure  of  a  composite  is  the 
maximum  stress  criteria.  The  failure  occurs  if: 

°L  —  °FN 

However ,  this  is  a  drastic  approximation,  since  the  fibers  vary 
significantly  in  their  strength.  Lee  proposes  that  in  addition  checKibg 
this  criterion,  the  fibers  fail  if 

where  OpS  is  the  fiber  shear  strength.  The  criterion  proposed  by  Hashin 
for  the  tensile  fiber  failure  is 

,2 
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Compression 

For  compressive  loads  applied  along  the  fiber  direction,  tbe  proposed 
failure  mechanism  is  analogous  to  the  buckling  of  a  column.  The  critical 
fiber  buckling  stress  in  the  shear  mode  is  given  by  Greszczuk  and  Hahn  as 

a  ■  G 
cs  r 

(1-k) 

where  Gr  is  the  resin  modulus,  and  k  is  the  fiber  volume  fraction  ratio. 
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Matrix  Failure: 

1.  Tension 

The  composite  tensile  strength  transverse  to  the  fibers  is  not 
expected  to  deviate  significantly  from  •'  e  matrix  tensile  strength. 
The  compressive  criterion  used  are  as  follows 

2  2  ** 

Lee  °t  -  °MN  or  (°TL  +  °TZ^  -  °MS 


TL  TZ'  -  ~MS 


Hashin  2  ^Z^  ^  ”  ^t^7^  ■>  Kt  ^  ®t7  ^  ^ 


where  (oT  *a^  )  >  0 
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2.  Compression 

Under  compression,  failure  may  occur  by  shearing  along  a  surface 
through  the  matrix  parallel  to  the  fiber  axis.  The  criterion  proposed  by 
Hashin  is 


°HHC  [(^s)  '  ‘l  (3TMZ)+ ^5—  (CtMZ>2  +  (0TZ“aT°Z)  + 
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where  is  the  compressive  matrix  strength. 
Delamination: 

Lee  proposes  that  delamination  occurs  it  either 

1  n  *S 


°Z  ±  °DN 


°r  (aLZ2+aTZ2)  >°DS 


where  and  o_g  are  the  through-the-thickness  tensile  and  shear  strength 
respectively,  xet  another  form  used  to  identity  delamination  is  as  follows 

/°z  \2  °lt2  +  ctTZ2  .  . 
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Damage  Prediction  Calculations 

The  damage  histories  tor  selected  coaposite  laminates  subjected 
to  both  in-plane  and  bending  loads  have  been  determined.  Mote  that 
the  loads  are  statically  applied.  Results  are  summarized  below. 
Uniaxial  Tension 

The  one  element  plate  model  of  Figure  25  is  employed  for  the 
uniaxial  tension  analysis  of  angle-ply  laminates.  The  laminate 
consisted  of  twenty-four  layers  of  T300/5208  graphite  epoxy.  The 
assumed  material  and  the  strength  properties  of  T300/5208  are 
given  in  Table  11.  The  first  and  the  last  ply  failure  curves  as  a 
function  of  lamina  orientation  angle  are  shown  in  Figures  26  and  27. 
As  expected  the  first  and  the  last  ply  failures  for  uniaxial 
laminates  of  [0/0/0]  of  Figure  2$  occur  simultaneously  where  as 
for  angles  greater  than  30,  they  are  quite  separated.  The  [0/90/0] 
layup  of  Figure  27  shows  that  for  0-60°,  75°,  and  90*  laminates, 
the  initial  and  final  failures  coincide  where  as  for  angles  less 
than  60°,  they  are  easily  distinguished.  Table  12  displays  an  alter-* 
nate  view  of  the  damage  progression  where  initial  failure  occurred 
at  the  second  load  increment  and  the  final  failure  at  the  fourth 
increment.  The  column  headings  of  Table  12;  TF,  CF,  TM,  CM  and  DL 
denote  Tensile- Fiber  Failure,  Compressive  Fiber  Failure,  Tensile 
Matrix  Failure,  Compressive  Matrix  Failure  and  Delamination 
respectively.  Thus  one  can  easily  identify  the  failure  mode 
within  a  ply  for  a  given  load  increment. 


Four-point  Bending 


The  bending  problem-  of  Figure  28  is  modelled  with  four  elements. 

The  material  and  strength  properties  are  as  listed  in  Table  11.  The 
laminate  is  unidirectional  and  consists  of  twenty  four  layers.  In 
the  bending  problem,  the  critical  aspect  ratio  is  defined  as  S*o 
Delamination  is  observed  f ot  aspect  ratios  less  than  the  critica 
middle  of  the  laminae  as  the  load  is  increased.  Additional  matrix  and 
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fiber  failure  accompany  delamination  as  shown  in  Table  13.  The  inter¬ 
action  curve  of  Figure  29  reveals  that  the  final  failure  occurs  after 
twenty  percent  load  increase  over  the  initial  failure  load.  It  should  be 
noted  that  for  aspect  ratios  less  than  the  critical,  the  percent 
increase  of  the  final  failure  to  that  of  Initial  failure  load  is 
constant;  thus  if  one  reduces  the  shear  strength  by  the  same  percentage, 
the  final  failure  load  for  the  corresponding  aspect  ratio  ends  up  on 
the  interaction  curve.  This  phenomenon  is  illustrated  by  the  dash-lines 
of  Figure  29.  When  the  aspect  ratios  are  higher  than  the  critical,  the 
fiber  failure  at  the  outermost,  laminae  proceeds  rapidly  toward  the 
center  and  within  four  percent  of  the  initial  load,  ultimate  laminate 
failure  occurs.  A  typical  damage  progression  is  displayed  in  Table  14 
for  an  aspect  ratio  of  100. 
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While  not  preeentea  herein,  some  work  has  been  performed  to 
relate  change  in  strain  energy  to  damage  progression.  It  is  hoped 
that  this  work  can  be  extended  to  calculating  strain  energy  release 
rate  as  a  function  of  the  extent  of  the  delaminated  region. 

11.3.3  TASK  111:  Correlation  of  Formulated  Response  Model  with 
Experimental  Data 

Some  quantitative  data  relating  to  the  impact  damage  of  composite 
specimens  has  been  assembled  [68-731.  It  will  be  utilized  along  with 
any  additional  data  obtained  to  perform  analysis/ test  correlations. 
Since  the  nonlinear  formulation  including  damage  effects  is  not 
complete,  the  only  use  of  test  data  has  been  of  that  in  [A3J  and 
discussed  in  the  previous  section. 
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Table  1.  Spurious  Zero  Energy  Modes  of  the  QHD  Family 
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TABLE  5 

NORMALIZED  STRESSES  AND  DISPLACEMENTS  FOR  SIMPLY  SUPPORTED  RECTANGULAR  0-90-0  PLATE  WITH  SINUSOIDAL  LOAD 
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TAuLE  6 

NORMALIZED  STRESSES  AND  DISPLACEMENTS  FOR  SIMPLY  SUPPORTED  SQUARE  0-90-90-0  PLATE  WITH  SINUSOIDAL  LOAD 
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Summary  of  Numerical  Examples 


Aspect 

Ratio 

Finite  Element  Solution 

Closed  Form 
Solution 
Reference  [4] 

QHD28 

QH040 

WBM 

2 

5.860 

5.525 

5.824 

5.500 

4 

9.780 

9.757 

9.706 

9.359 

10 

15.440 

15.340 

15.276 

15.145 

20 

17.850 

17.719 

17.628 

17.665 

25 

18.246 

18.103 

18.006 

18.093 

100 


18.964 


18.805 


18.704 


18.733 


TABLE  9 

The  Effects  of  Linear  (z),  Quadratic  (z2)  and  Cubic  (z3) 


Table  11  Material  Properties 


Elastic  Constants 

Uniaxial  Tension  Four-point  Bending 


Et  (GPa) 

138 

190 

E2  (GPa) 

10.6 

11 

E3  (GPa) 

10.6 

11 

G12  (GPa) 

6.4 

7.2 

Gu  (GPa) 

6.4 

7.2 

G23  (GPa) 

6.4 

7.2 

V12 

0.3 

.38 

V13 

0.3 

.38 

V23 

0.3 

.38 

Strengths 

0FN  (MPa) 

1500  (1500)* 

1502  (1502) 

ops  (MPa) 

68 

67.5 

°MN  <»•> 

40  (246) 

41' (250) 

°MS  (MPa) 

68 

67.5 

°DN  <»«> 

40 

41 

aDS  <MPa> 

68 

67.5 

♦Terms  in  parenthesis  are  the  compressive  strength. 


Table  12  Damage  Accumulation  of  a  24-ply, 
[30/90/-30]  laminate  under 
uniaxial  load. 


M.T  If  Cf  TM  CM  01 

1  0  0  2  0  0 

2  0  0  2  0  0 

3  0  0  2  0  0 

4  0  0  2  0  0 

5  4  0  0  0  0 

6  4  0  0  0  0 

7  4  0  0  0  0 

8  4  0  0  0  0 

9  0  0  2  0  0 

10  0  0  2  0  0 

11  0  0  2  0  0 

12  0  0  2  0  0 

13  0  0  2  0  0 

14  0  0  2  0  0 

15  0  0  2  0  0 

16  0  0  2  0  0 

17  4  0  0  0  0 

18  4  0  0  0  0 

19  4  0  0  0  0 

20  4  0  0  0  0 

21  0  0  2  0  0 

22  0  0  2  0  0 

23  0  0  2  0  0 

24  0  0  2  0  0 
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Damage  Progression  of  a  24-ply  laminate 
with  Aspect  Ratio  S-8  under  the  four- 
point  bending  load.  W1-20- 
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Figure  1:  SPURIOUS  ZERO  ENERGY  MODES  (QHD23  ELEMENT) 
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NORMALIZED  TRANSVERSE  DISPLACEMENT  VS.  ASPECT 
10  FOR  CYLINDRICAL  BENDING  OF  BIDIRECTIONAL  AND 
SYMMETRIC  3-PLY  LAMINATES 


NORMAL 


NORMALIZED  IN-PLANE  STRESS  FOR  C 
OF  SYMMETRIC  LAMINAT 


FEM(  PRESENT) 


ZEO  SHORT-TRANSVERSE  NORMAL  STRESS  FOR 
PPORTEO  0-90-0  SQUARE  PLATE  x=a/2,  yaQ 


NON-DIMENSIONAL  I  ZED  FUNDAMENTAL  FREQUENCY  VS.  ANGLE 
OF  PLY  ORIENTATION  FOR  ANGLE-PLY,  SIMPLY  SUPPORTED, 
SQUARE  PLATE  OF  PROBLEM  3 
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FIGURE  19: 


NON-DIMENSIONALIZED  FUNDAMENTAL  FREQUENCY  VS.  ANGLE 
OF  PLY  ORIENTATION  FOR  ANGLE-PLY,  CANTILEVER  PLATES 

OF  PROBLEM  4  ,  , _ 
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DISPLACEMENTS(CM‘  10D-3) 


tint  (HlCROSteOMDS) 

FIGURE  21:  ONE-MODE  DISPLACEMENT  VS.  TIME  RESPONSE  OF  A  COARSE-MESH 

(0/90)  LAYUP  SQUARE  PLATE  UNDER  SUDDENLY-APPLIED  SINUSOIDAL 
LOADING 


«ure 27.  Fir3t  and  Last  Ply  Failure  Curves  for  (0<7S9#/-e]  laminates 
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APPENDIX  IA  -  HIGHER  ORDER  DISPLACEMENT  MODELS 


QHD40 


NODAL  DEGREES  OF  FREEDOM: 


Corner  Nodes  -  (u0  v0  w0  ity  4>x  <J)y}T 
Mid-side  Nodes  -  {w0  \px  ipy}  7 


DISPLACEMENT  FIELD: 


U  =  U0  +  Z\px  +  Z2<t>x 


V  =  V0  +  Z^y  +  Z  2<t>y 


W  =  Wr 


where; 


T 

uo»  v0»  ^X»  4>y  :  U  X  y  xy}  {a} 

w0»  ’>x»  'Py  :  U  x  y  xz  xy  y2  x3y  xy3}7  {B} 


STRESS  FIELD: 


From  constitutive  relations  -  0|  *  C-jje-jj  (orthotropic  mat.) 
axx  =  f(z2,  x2,  y2) 

°yy  =  f(z2»  *2>  y2) 
oXy  s  f(z2»  x2,  y2) 

From  equilibrium  considerations  -  a-,* j , j  =0 
<?xz  s  f(z3.  x  *  X ) 

0y2  s  f(z3,  x,  y) 

<7ZZ  *  f(z3) 
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NODAL  DEGREES  OF  FREEDOM: 


QHD28 


*uo  vo  wo  Vy 


DISPLACEMENT  FIELD: 


U  +  U0  +  2I|JX  +  Z2<tx 
V  =  VQ  +  ZiPy  +  Z2<t>y 


where ; 


w  =  w0 

u0‘  vo'  V  '^x>  ^y *  <J>X*  <f>y*  :  u  X  y  xy}T{a} 


STRESS  FIELD: 

i.  From  constitutive  relations  -  o-,-  =  C^jeij  (orthotropic  mat 
°xx  =  ^(z2,  x,  y) 

°yy  =  f(z2»  X,  y) 

°xy  =  f(z2.  x,  y) 

11.  From  constitutive  considerations  -  aij,j  =  0 

0X2  =  f(23) 

OyZ  *  f(z3) 

azz  =  constant 
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APPENOIX  IB  -  MODIFIED  KIRCHHOFF  FORMULATION 

QD32 

NODAL  DEGREES  OF  FREEDOM: 

Corner  Nodes  (w0  v0  w  3—  yx  Yy} 

Mid-Side  Nodes  (w  |H}T 

dn 

DISPLACEMENT  FIELD: 

w  =  f(x,  y) 
u  =  uo  “  z(|^  +  Yx) 
v  *  vo  -  z(l7  +  Yy) 

where ; 

uo*  vo>  x1  y  :  {1  x  y  *y}T{a} 

w  =  { I  x  y  x2  xy  y2  x3  x2y  xy2  y3  x1*  x3y  xy3  y4  x“y  xy4}T{B> 

STRESS  FIELD: 

i.  From  constitutive  relations  -  o ^  =  CiJ£ij  (orthotropic  mat.) 
oxx  =  f(z,  x2,  y2) 

<Jyy  -  f(z,  x2,  y2) 
aXy  =  f(z,  x2,  y2) 

ii.  From  equilibrium  considerations  -  o^j  =  0 

°xz  =  f(z2»  *2*  y2) 
oy2  =  f(z2,  x2,  y2) 

ozz  =  f(z3,  x,  y) 


w 

w|2^YxY„)T  ^ 


NODAL  DEGREES  OF  FREEDOM 


{uO  VoW^^YxYy) 

(w  l^}T 

cH 


Corner  Nodes 
Mid-side  Nodes 
Center  Node 


DISPLACEMENT  FIELD 


w  =  f(x,  y) 


C  ?W  <3W  V  Y 

tun  vn  W  T-  —  x  ' 
0  0  ?X  - y 

<„>T 

{w  ^  lw>T 

dX  dy 


(I1:  *  **) 


where; 


v  =  vo  '  2 


V  V  x,  y  :  U  x  y}T(a} 


w  :  (1  x  y  x2  xy  y2  x3  x2y  xy2  y3  x“  x3y  x2y2  xy3  y*}‘ 
STRESS  FIELD 

i.  From  constitutive  relations  -  cr^  =  C^e^  (orthotropic  mat.) 
axx  =  f(z,  x2 ,  y2 ) 

Vyy  =  f(z,  x2 ,  y2 ) 
tJxy  =  f(z-  x2 ,  y:  ) 

ii.  From  equilibrium  considerations  -  0 i j » i j  =  0 
axx  =  f(z2,  x,  y) 
ayz  =  f(z2,  x,  y) 
ozz  =  f(z3) 


V  S  \  A  *. 


APPENDIX  II 


MASS  MATRIX  FORMULATION 


The  mass  matrix  for  elements  under  development  is  easily  arrived  at 
by  considering  kinetic  energy  in  the  form 

T  *  j  ^p(u2  +  v*  +  w2)dV 

V 

where  u,  v  and  w  represent  displacements,  p  is  the  mass  density  and  the 
dot  superscript  denotes  velocity.  Defining  velocities  in  terms  of  element 
shape  functions  gives 

T  «  \  Ca]T  J.P(NU}[NU]  +  (NV}[NV]  +  {NW}[NW]  dV  {A> 

V 

which  is  the  classical  form 
T  =  \  [a]T[M]{a) 

The  element  mass  matrix  [M]  is,  therefore,  specified  as 

CM]  =  J  p{Nu}{Nu}  +  {NV}{NV}  +  {NW}{NW}  dV 

V 

Note  that  the  shape  functions  [N^]  involve  distance  from  the  mid-plane 
of  the  element  to  a  layer  denoted  by  Z  and,  therefore,  the  mass  matrix 
definition  provided  not  only  represents  mid-plane  inertial  effects  but 
also  rotatory  inertia  as  well. 


APPENDIX  III  -  LARGE  DISPLACEMENT  FORMULATION 

Based  on  Green's  Strain  Tensor,  the  following  procedure  is  utilized 
to  obtain  the  large  displacement  and  the  geometric  stiffness  matrices. 

Let  N  be  shape  functions  relating  displacements  at  any  point  in  the 
element  {6}  to  nodal  displacements  {A}  such  that 

{5}  =  [N]{a> 


Also  let  (N-j,j}T  denote  those  shape  functions  associated  with  the  ith 
displacement  field  (i  -  u,v,w)  and  ",j"  denotes  the  differentiation  with 

respect  to  the  jth  coordinate,  i.e.,  where  x,  =  x,  x2  =  y  and  x3  =  z. 

°x  • 

Then,  the  strain  exx  given  by  J 


=  ly.  + 1  /Iy.Y+  /  Y+  ( ?wY 

3x  2  \3 x/  \3x  /  \3x/ 


can  be  written  as 


=  ^U’X^  *•-}  +  2  l^U’X-  ~^U’X;  +  1-^V»X;  ^V’X^  +  ^W>X;  ^W»x-'l  \ii} 


Similarly;  the  shear  strain  eY,.  can  be  represented  by 

eXy  =  {{Nu,y}T  +  {Nv,x}T]-{A}T  +  {A)T({Nu,x}T{Nu,y}  +  {Nv,x}T{Nv,y} 


+  {Nw,x}‘{Nw,y}}  {A} 


The  strain  field  in  indicia!  notation  is  expressed  by 

{eij}  =  l|[{Ni’j}T  +  {NJ*i}T]{A}  +  (A>T[tNk»i>T{Nk,j}]{A} 
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Then  the  incremental  representation  becomes 


{5eij}  *  l[{{Ni’j}T  +  {MJ»i}T}{t’A}  +  t6^T[{Nk.i>T(Nk,j}]{A} 
+  {A}T  [{N|t,1-}T{Nk,j}]  {6A}j 


But  the  second  term  can  be  expressed  as 


{A}T[{Nk,j}T{Nk,i}]  (SA) 


Thus  combining  terms 


CB0]  =|[{Ni‘j}T+  {NJ‘i}T] 

[Bl]  =  ±{A}T  [{Nk,i}TlNk,j}  +  {Nk>j}T  {Nk,,}]  = 


{A}|[Mxx] 

(A}T[Myy] 

{A)^[MXy] 

(A}T[Mzz] 

{a}t[mxz] 

(A}T[Myz] 


.<v,v 


. *  -  ■  ** M r 


Uelj)  *  UftNl.j)7  +  (Nj.1>T](«M  ♦  UC  [(Nkl1}T{Nk,j}  +{Nk,j)T{Nk„l]  {iA}J 

fcS&s 


■ 


i*r^ 

ft.. 


(6A^j)  *  [B0]{6A>  +  [BL]{<SA} 


where  [B0]  is  the  linear  component  and  [Bl]  is  the  large  displacement  component 
Having  the  definitions  for  [BQ]  and  [Bl3»  the  small  and  large  displacement 
matrices  [Kq]  and  [KL]  are  represented  as 

[K0]  *  J  [B0]T[D]CB0]dV 

V 

W  *  J  |cbl]t[0][bo]  +  [BL3T[0][BL]  +  [Bo]T[0][BL]JdV 

V 


The  geometric  stiffness  matrix  is  also  derived  from  [BL]  and  it  has  the 
following  form 


[KG]  *  J  (oxxCMxxJ  +  OyyCMyy]  +  czz[ Mzz]  +  oxy[Mxy]  +  oxz[MX2] 


Where  the  o's  are  the  stress  components  and  again  integration  is  on  a 
layer  by  layer  basis. 


